First-principles investigation of dynamical properties of molecular devices under a steplike pulse 
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We report a computationally tractable approach to first principles investigation of time-dependent current of 
molecular devices under a step-like pulse. For molecular devices, all the resonant states below Fermi level 
contribute to the time-dependent current. Hence calculation beyond wideband limit must be carried out for a 
quantitative analysis of transient dynamics of molecules devices. Based on the exact non-equilibrium Green's 
function (NEGF) formalism of calculating the transient current in Ref@, we develop two approximate schemes 
going beyond the wideband limit, they are all suitable for first principles calculation using the NEGF combined 
with density functional theory. Benchmark test has been done by comparing with the exact solution of a single 
level quantum dot system. Good agreement has been reached for two approximate schemes. As an application, 
we calculate the transient current using the first approximated formula with opposite voltage Vt(t) = -Vu(f) 
in two molecular structures: AI-C5-AI and Al-Cf,o-Al. As illustrated in these examples, our formalism can be 
easily implemented for real molecular devices. Importantly, our new formula has captured the essential physics 
of dynamical properties of molecular devices and gives the correct steady state current at t = and t — > 00. 

PACS numbers: 71.15.Mb, 72.30.+q, 85.35.-p 73.23.-b 



INTRODUCTION 

With the rapid progress in molecular electronics, [Q] quan- 
tum transport in molecular device has received increasing 
attention. In particular, the dynamic response of molecu- 
lar devices to external parametersJll-lgt], in which the exter- 
nal time-dependent fields or internal parametric pump poten- 
tials drive the electrons to tunnel through the molecular de- 
vice, is one of the most important issues in molecular elec- 
tronics. The simplest molecular device structure is the two- 
probe lead-device-lead (LDL) configuration, where "device" 
is the molecular device connected to the external probes by 
the "leads". In such a device, all the atomic details of the 
device material can be treated using density functional the- 
ory (DFT) and the non-equilibrium physics can be taken into 
account using non-equilibrium Green's function (NEGF). Up 
to now, from an atom point of view, one of the most popu- 
lar theoretical approaches used to study the quantum trans- 
port properties of molecular device is Keldysh nonequilib- 
rium Green's functions coupled with density-functional the- 
ory (NEGF-DFT).|li Using this approach, the steady state 
quantum transport properties in molecular devices have been 
widely studied. 

For time dependent response of molecular devices, there 
have been many different theoretical approaches, such 



as evolution of time-dependent Schro ding er equation, yjj, 
time development operator approach, [12] and the NEGF 
technique. ll 1311 These approaches are convenient to deal with 
dynamic response of time-dependent external field that is si- 
nusoidal (e.g., microwave radiation). Under such an external 
field an electron can tunnel through the system by emitting 
or absorbing photons, giving rise to the photon-assisted tun- 
neling (PAT). Concerning the steady state ac response to har- 
monic external field, the Floquet approach is convenient. 114 1 
For the transient transport, however, the pulse like ac sig- 



nal is the optimal driven force since they can provide a less 
ambiguous measure of time scales. 11511 In this case, besides 
PAT, one of the most interesting questions to ask is how fast 
a device can turn on or turn off a current. With the de- 
velopment of molecular electronics, providing a particular 
viable switching device has become a key technical issue. 
Concerning the transient dynamics, different approaches such 
as path-integral techniques, (3 the solution of Wigner dis- 
tribution function,|[n| the time-dependent numerical renor- 
malization group, jlif time-dependent DFT (TDDFT),|3l[l9ll 
and Keldish Green's funcitonfl 0, [20I] have also been de- 
veloped and applied to different systems. Up to now, most 
of these approaches can only be implemented in simple sys- 
tems such as quantum dotsj^l [2(J or one-dimension tight- 
binding chains. Numerical calculation of transient current 
for molecular devices is very difficult at present stage due to 
the huge computational cost. This is because if we calculate 
the current as a function of time t, the amount of calculation 
scales as f 3 if the time-evolution method is used. This scaling 
can be reduced to a linear scaling in t\f the wideband limit is 
used. [ 21 ] As we have demonstrated,^! the wideband limit is 
not a good approximation for molecular devices. If one uses 
the exact solution from NEGF,|0] one can calculate the tran- 
sient current at a particular time. However, the calculation in- 
volves a triple integration over energy which is extremely time 
consuming. Clearly an approximate scheme that is suitable for 
numerical calculation of transient properties for real molecu- 
lar devices while still captured essential physics is needed. 

It is the purpose of this paper to provide such a practical 
scheme. To study transient dynamics, in this paper, we con- 
sider a system that consists of a scattering region coupled to 
two leads with the external time dependent pulse bias poten- 
tial V a (t) = 6(+t)V a . For this case, the time-dependent cur- 
rent for a step-like pulse has been derived exactly without us- 
ing the wide-band limit by Maciejko et al|0]. Since the gen- 
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eral expression for the current involves triple integrations, it 
is extremely difficult to perform them in a real systems like 
molecules devices. So, approximation has to be made. The 
simplest approximation is the so called wide-band approxima- 
tion where self-energies U- a are assumed to be constants inde- 
pendent of energy. 12311 Unfortunately, this approximation can 
not give the correct result since in general there are several res- 
onant levels that significantly contribute to the transient cur- 
rent in molecules devices. To go beyond the wideband limit, 
we propose an approximate scheme of calculating the tran- 
sient current that is suitable for numerical calculation of real 



molecules devices. 12411 Our scheme is an approximation of the 
exact solution of Maciejko et al[6]. It is very fast computa- 
tionally and gives the correct limits at t — and t — oo. Since 
the exact solution of transient current is available for a single 
level quantum dot system, we have compared our result with 
the exact solution on the quantum dot system to test our ap- 
proximate schemes. Good agreement is obtained. Therefore, 
our approximated scheme maintains essential physics of tran- 
sient dynamics. Using our scheme, we calculate the transient 
current for the upward pulse (turn-on) and downward pulse 
(turn-off) in two molecular structures: AI-C5-AI and AI-C60- 
Al. We find that different from the single level quantum dot 
system, upon switching on the current oscillates rapidly in the 
first a few or tens fs with several characteristic time scales. 
Furthermore, due to the resonant states in molecular devices, 
transient currents have a much longer decay time r, especially 
for the molecule device having a complex electronic structure 
such as AI-C60-AI. 

The rest of paper is organized as follows: In Sec. II, starting 
from the typical molecular device Hamiltonian which is ex- 
pressed in an non-orthogonal basis, we shall derive a general 
DC and AC current expressions for a non-orthogonal basis 
set. It is found that for DC bias, the expressions of current 
for orthogonal and non-orthogonal basis sets are the same. 
For ac current, however, the expressions are different as will 
be demonstrated in Sec. II. The reason that we study the dif- 
ference between orthogonal and non-orthogonal basis sets is 
the following. For the NEGF formalism, it is assumed that 
the basis set is orthogonal. It turns out that for ac transport, 
the current expression becomes extremely complicated if non- 
orthogonal basis is used. For DFT calculation, however, most 
people work in molecular orbitals that are non-orthogonal. 
Our results show that we must orthogonalizing the nonorthog- 
onal molecular Hamiltonian, so that the present approach in 
RefH can be used. In Sec. Ill, based on the exact solution of 
Maciejko et al, we derive two approximate expressions for 
transient current with different levels of approximation. They 
are all suitable for numerical calculation for real molecular 
devices. In addition, the initial current and its asymptotic long 
time limit are shown to be correct. In Sec. IV, in order to appre- 
ciate our approximate formulas, we compare our result with 
the exact result obtained in Ref.6 for a single-level quantum 
dot connected to external leads with a Lorentzian linewidth. In 
Sec.V, we apply our formalism to several molecular devices. 
Finally, a conclusion is presented in Sec. VI. Two appendices 



are given at the end of the paper. In Appendix A, we give a 
detailed derivation of orthogonalization relation for an non- 
orthogonal basis. This relation is used to derive the effective 
Green's function which is the key to approximate exact cur- 
rent expression of Maciejko et al. In Appendix B, we show 
how to orthogonalize an nonorthogonal Hamiltonian so that 
the general AC current for real molecules device can be de- 
rived. 



GENERAL AC CURRENT 
Hamiltonian 

The transport properties of a molecular device can be de- 
scribed by the following general Hamiltonian: 



H = H C + H T + J] H a 



(1) 



a=L,R 



where Hi and Hr describe the left and right macroscopic 
reservoir, respectively; H c is Hamiltonian of the central 
molecular device; Hj couples the reservoirs to the molecu- 
lar device. For a particular basis set, the above Hamiltonian 
can be written in the following matrix form: 

v a ,v c 

where e is the electron charge, c Va (c y t) and d Vc (d\ c ) are 
Fermionic annihilation (creation) operators at the state v in 
the lead-a and the state v in central molecular device. v a , 
v c are the indices of the given basis set. The Hamiltonian of 
lead-a are divided into two parts: the time independent part 
and time dependent part due to external bias V a (t) on the 
lead-a. Here we consider two kinds of step-like bias: up- 
wards pulse (turn-on case) V„(f) and downwards pulse (turn- 
off case) V® (t), where 



„ ( , V a , t<0 
a K > 10, t > 



In the adiabatic approximation it is assumed that the sin- 
gle particle energies acquire a rigid time-dependent shift as 
H° + W a (t). The energy shift in the leads is assumed to be 
uniform throughout. This assumption is reasonable since the 
pulse rising time is slower than the usual metallic plasma os- 
cillation time, which ensures that the external electric field is 



effectively screened. 12511 

Since Green's function G r (t, /') is obtained by solving 
Dyson equation from the known history, it is better to set 
time dependent external bias V a (t > 0) = so that the un- 
certainty of future can be eliminated. Ja] From Eq.©, this 



3 



is satisfied only in the downward case. In the following, 
we will discuss how to eliminate this uncertainty for the up- 
ward pulse. To use the Dyson equation, we will separate the 
Hamiltonian into two pieces: the unperturbed Hamiltonian 
that can be exactly resolved and the interacting term which 
contributes to the self energy in Dyson equations. For the 
downward pulse, we define the non-biased open system as 
the unperturbed system. It is described by the Hamiltonian 
H° = H° + Hj! + Hj. For the upward pulse, however, the sit- 
uation is different, in which we will set the DC biased open 
system H v = [H° + V a l] + [H° + U v ] + as the unperturbed 
Hamiltonian and set V^(t) = V^(t) - V a as the new time de- 
pendent part. Here denotes the coupling between scatter- 
ing region and biased leads and U v is the induced coulomb 
potential due to the external bias. Now, the time dependent 
bias V u satisfies V u (t > 0) = 0, and the uncertainty of the fu- 
ture in the upward case is eliminated. Then, for the downward 
case, we have V°(t) = V°(i) and W = H° while for the up- 
ward case we have V%(t) = V%(t) - V a and H" = H v . From 
now on we will use superscript "ex" to denote the unperturbed 
system that is exactly resolvable. 

When the system is biased, the incoming electron will po- 
larize the system. The induced Coulomb potential in the cen- 
tral scattering region consists of two parts: DC and AC parts. 
The DC part can be put into the exactly resolvable Hamil- 
tonian H". The induced time dependent coulomb potential 
U(f) due to the external bias V a (t) is included as part of the 
non-equilibrium Hamiltonian. Because the electric field is 
not screened in the small scattering region where the poten- 
tial drop occurs, the coulomb potential landscape U(t) in the 
central region is not uniform, which is different from the semi- 
infinite leads. Note that it is rather difficult to treat the time- 
dependent coulomb potential and no close formed solution 
exists if one does not assume wide band limit. In the small 
bias limit, we can expand the time-dependent coulomb po- 
tential to linear order in bias U(f) = e2 Q ,u Q ,V' Q .(f) so that the 
analytic expression for current can be obtained. Here u a is 
the characteristic potential. 12711 From the gauge invariance, 
ll26ll Xn-Up. = I, and u„ is determined from a poisson like 
equation. ||28| In this paper, we consider the symmetric cou- 
pling so that for the external bias Vz.(f) = -V«(f) it is a good 
approximation to assume that the time dependent coulomb po- 
tential U (f) is roughly zero in the the molecular device regime. 

In the following, we will derive an exact solution of tran- 
sient current using a non-orthogonal basis set. J29h To facili- 
tate the derivation, we take a unitary transformation (9(f) to 
the Hamiltonian (O with 



0(t) = expiieJ^j^dT [v a (T)cl a c Ya ] 



where V a (r) = 6(-T)V a for the downward pulse and V a (T) = 
—6(—T)V a for the upward pulse. Note that the time f in (9(f) 
can be negative or positive, and (9(f) = 1 only when f > 0. The 
new Hamiltonian 'H = 6Hd\t) + i{^.d{t))&{f), in which "H a 



and Irij are different from original ones: 

a / I jig ftgVg V a 

•Hr = J^clT VaVc (t)d Vc + h.c. 



(4) 



where 



T VnV( .(f) = TS BVe SB„(/) 

3B„(?) = expfj'e I V a (r)dT] 
Jo 



(5) 



For the original Hamiltonian with nonorthogonal basis, the 
overlap between nonorthogonal basis is expressed as the ma- 
trix form Sj] v = (ji\v). After the unitary transform, annihilation 
(creation) operators c a (c„) and consequently the orbital basis 
/i„ in the leads are changed, then overlap matrices between the 
leads and the scattering region become 



S VaVc (t) = s° VaVc m a (t) 
s YeVa (t) = ml(f)S° VcVa . 



(6) 



In the following, we will use the transformed Hamiltonian 
[Eq. ( 1415b . in which c Va , d Vc are used] to derive the time de- 
pendent current expression. 



The current 

The current operator from a particular lead-cr to the molec- 
ular junction can be calculated from the evolution of the num- 
ber operator of the electron in the semi-infinite lead-o-. As- 
suming there is no direct coupling between the left and right 
leads, the current operator can be expressed as: lf30ll 



lit) = -eYj- f ^ v " (t) 

= -'Z 

= e J] c\ a {t) \iT VaVc (t) + S Vc (f)~ J d Yc (t) + H.c.(l) 

where 'He' denotes the Hermitian conjugate. The current is 
obtained by taking average over the nonequilibrium quantum 
state '< ... >', 



J a (t) = e J] 



T Va , Vc (t')-S Va , Vc (t')i-\ 



(8) 



/=/' 
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where "J^" and "Jj" denotes the left and right derivation re- 
spectively, and 

G< aVa (f, t') = i (cl(t')d Vc (t)) , G< nVt {tf, t) = i «.(f)c Vn (f')) ■ 

Using the Keldysh equation and the theorem of analytic 
continuation, we have 

G< a (t,t') = fdti [G r cc (t,t l )B ca (t l )g < aa (tut')+ 

G f < c (f,f 1 )B Cff (f 1 )g^(f 1 ,f')] (9) 



where 



B ra (fi) - T ra (fi) - S ra (fi)/— 

at 



(10) 



For simplicity, we have dropped the subscript ji, and keep only 
the symbol c and a to indicate the central scattering region and 
lead-o-, respectively. In the above expression and in the fol- 
lowing, the summation convention on repeated sub-indices is 
assumed. Substituting Eq.© into Eq.©, we have the general 
expression for the current: 

J a (t) = -2eRe j dh Tr 

[G r cc (t,h)B ca (h)g2 a (h,f)Kc(t')- 

G%(f, h)B ca (t 1 )g a aa (t l ,f)B ac (t')] t=t , (11) 

When the system reaches a stationary state, V a (f) - V a be- 
comes time independent, from definition Eq.Q, and ( fTUl ), 
we can find 

n ca (h)XB ac (t) = e- ieV ^-^ ca ml c , 

with B° , = T° . - i4-S° , , where "0" denotes the zero 

ca/ac ca/ac at ca/ac 

bias system.In addition, all the propagators G and g depend 
only on the time difference t\ - t. Taking the Fourier trans- 
formation, from Eq.® or Eq. dTTb . we can easily obtain DC 
current expressed in the energy representation: 

J a = J dej a {e) 

= Re2eJ de Tr [G r (e)E<(e) + G^E^e)] (12) 

where G and T, are the Green's function and the self-energy. 
They have the same matrix dimension as that of the Hamil- 
tonian H c . The Green's function G^" and self-energy Yf^ a is 
defined as 

G rla {e) = [eI-H c -£' /fl (e)] -1 

11(e) = [T° a - e a S° ca ] gl a (e a ) [T° c - e a Sl ] ( 1 3) 

where e" — e - eV a , I is the unitary matrix with same dimen- 
sion as H c , y — r,a, <, and 



((e±iO + )S aa -B° aa )' 



is the surface Green's function of the semi-infinite periodic 
lead which can be calculated numerically using a transfer ma- 
trix method. 13 ill Here, f(e) is the Fermi distribution. Eq. (fT2b 
shows that the dc current expressions are the same for both 
orthogonal and non-orthogonal basis sets. 

When the time dependent field V„(t) is present, however, 
the current expressed in energy representation will be very 
complicated for nonorthogonal basis due to the term S(t')ij t 
in Eq.®, since B(?i)XB(f) can't be expressed as a function of 
time difference t\ - t. One thing is clear, the transient current 
expressions are different for orthogonal and non-orthogonal 
basis sets. Instead of deriving a complicated transient current 
expression using a non-orthogonal basis set, we will elimi- 
nate S ca / ac (t')i^ in Eq.® and work on an orthogonal basis 
set. In Appendix D from the overlap matrix S, we derive the 
orthogonal basis set and new Hamiltonian H expressed in this 
orthogonal basis. With the new orthogonal Hamiltonian, the 
overlap matrix S ca / ac (f) will be eliminated since the overlap 
matrix of orthogonal basis S°"'' = I. Then, replacing Hamilto- 
nian H in Eq.© with H and go through the derivation leading 
to Eqs.(t2lfTTT> again, we arrive at a new AC current expression: 



/„(*) = 2eRe J" dhTr {G r cc (t, h) [T ca (ti)gtf(h - t)T ac (t)]\ 
+2eRe j dt 1 Tr{G< c (t, h) [T ca {h)g a a e a x {h - t)T ac (t)]} (15) 

Defining the self-energy on the orthogonal basis 



^y=r,a,< 



(t,t') = T ca {t)gZ x {t-t')T ac {t') 



(16) 



V a; €SUr,jU a .£SlH* 



g^r(e) = f{e)[g a aa {e)-g' aa {£)} 



(14) 



where gaa\t-f) = f || e^' c( -'^'" > g y ^ x {e) is the surface Green's 
function of semi-infinite lead-ff in the unperturbed state as de- 
fined in the Sec„ For the downward pulse we have set the 
unperturbed system as the open system at zero bias, in which 

gaa X (f ) = U - + ;0 + . For the upward pulse, the un- 

L Jo-esur 

perturbed system means V a biased open system, in which 

gla\e) = U- eV a - H l l + /0 + l 1 . From Eq.fB),©, we 
have the general current formula 



/„(/) = 2eRe J" dt<Ir [G r (t, h)^(h, + G K (t, h)I^{tu 0] 

(17) 



At t < 0, AC external bias V a (t) or time dependent part in 
Hamiltonian V a (t) is a constant and the system is in a steady 
state. Consequently, the total current is known from DC trans- 
port theory that is expressed in the form of Eq.([T2b but with 
the Green's function and self-energy obtained from the or- 
thogonal Hamiltonian defined above. Hence in the following 
we shall derive only the Ac current when t > 0. First, we shall 
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look at the self-energy. From Eq.© and ( TTol l. 

rS(t,0 = ^w[T° ag L(*y)ry2B a (*') 



J ^ 



= 9Bt (oaJ t (f) | J* 



de 
2^ 



® a {t'W a {t') 
(18) 



where 23„(f) = 1 for the downward pulse and 23„(f) = e «' for 

y 

the upward pulse. Here H„ (e) is the self-energy at zero bias, 
- T° ,gJ ( "(e)T" [ . is the self-energy at the unperturbed 

In the 

7,0 



state defined above. In the downward case £„ 



upward case Hj' 



Setting S° ac = S° ca = 0, E£° and 
are defined in Eq.dlJil with zero and nonzero V a , respec- 
tively. We have Tl'^ie) = E^ /a '°(e - eV a ). From Eq.([T7]) and 
([T"8l , we find 

/„(*) = 2eRe J g £ dt x 

[G''(f, fi)^(e, ti.O + G < (», fi)^(e, fi, t)] (19) 

where the first term is the current flowing into the molecular 
device while the second one is the current flowing from the 
molecular device, and 



(20) 



y 

where 'Wait) = ^ a (t)W a (t). Here "L„ a is the self-energy of 
lead-a at zero bias. The lesser Green's function is given by 



G < (f, t 



■>-/*/ 



dh G r (t, t\) 



2>|(f 1( / 2 ) 



G a fe,f') 



f dt 2 e- ie «-»<W fs (t 2 )G a (f,t 2 )<Wl(t) 

%J — oo 



(21) 



Substitute Eq.d20b and (f2TT > into Eq.(TT9b and introducing a 
spectrum function 

A a (t,e)= f dt i e i «-V<W a {t)G r {t,t i )<Wl{h) (22) 

%J — oo 



we have 



J%(f) = 2eRe j g A„(t, (23) 
2eRe J ^ ^fce^^fU) (24) 



•C(0 



where 



(25) 



Very often, J7^ a (t - t') is singular at t — t', such as 
the quantum dot system with the wide-band limit S r ' fl (0) 
j" ^£ Ul a (E) = 5(0)(+r/2), or the superconducting-quantum 
dot-normal metal system, and so on. In these cases, we should 
be careful with Eq.d25t. 



Fg a (t,e) = F«„(f, e) + F« a (f, e) 



£4D— If 

Akf',e)nV^(f')K'°(£)^a) 



lE(f-r') 



(26) 



The first integral J ^ is the same as Eq.(l25b. the second inte- 
gral j J r now becomes Fp a (t, e) = A^(f, e)A°, where we have 



defined 



= i f e;,'" "(0) 

Then, Eq.(l24b becomes 

■£"*(*) = 2eRe J ^ ^A^t, f )i; ( e )FMU) 



(27) 



+ 2eRe 



/ £ £A^t,e)E;- (e)A*(t,e)A* (28) 



We note that Eq.(l28l> is the same as that derived in Refjfjl 
Different from Ref|g, we have split the expression into two 
terms. The first term corresponds to the non-wideband limit, 
i.e., when the linewidth function Y goes to zero at large en- 
ergy. The second term of Eq.d28T> is related to the wideband 
limit. Hence, for a quantum dot with a Lorentzian linewidth 
function^, only the first term is nonzero while for the system 
in contact with a superconducting lead both terms are nonzero. 

So far, we have discussed the ac conduction current J a (t) 
under the time dependent bias derived from the evolution of 
the number operator of the electron in the semi-infinite lead-a. 
Now we wish to address the issue of charge accumulation in 
the scattering region. In principle, this can be done by includ- 
ing the self-consistent Coulomb potential due to ac bias. l28ll 
However, at finite voltages, there is no close form expression 
for ac current if Coulomb potential is included. Alternatively, 
one can treat Coulomb potential phenomenologically as fol- 
lows. From the continuity equation, £ ff Ja(t) + dQ(t)/dt = 0, 
we see that the conduction current is not a conserved quan- 
tity. In the presence of ac bias, the displacement current J* 
due to the charge pileup dQ/dt inside the scattering region 
becomes important and must be considered. Since we have 
neglected the Coulomb interaction in our calculation, we can 
use the method of current partition H32l 13311 to include the dis- 
placement current. This can be done by partitioning the total 
displacement current Yia-fa - dQ/dt into each leads giving 
rise to a conserving total current I a = J a + J*. For symmet- 
ric systems like what we shall study below, it is reasonable to 
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assume that J d L = 7^ from which we find jf t = -{Jl + -/r)/2. 
Hence the total current is given by Ii = [Ji - /«)/2[25] which 
satisfies the current conservation I L + I R - 0. 



TRANSIENT AC CURRENT 

Up to now, we have derived the general expression for time 
dependent current, Eq. d22l23|25|28| > which can be used for or- 
thogonal as well as nonorthogonal basis set. To calculate the 
transient current we have to solve the retarded Green's func- 
tion G'(f, t') and integrate it over time to find A/?(/, e) and 
Fp a (t, e). For the pulse-like voltage V a (t) = +6(-t), we can 
obtain the Green' function G'(f, f) by solving Dyson equa- 
tion G r = G r - eq + G'~' eq EG r from the known history in the 
time domain. Depending on what is the chosen unperturbed 
system that can be solved exactly, the Dyson equation can be 
written in a different but equivalent form. In the study of time- 
dependent transport, it is better to treat the time-independent, 
open steady state system as the unperturbed system as de- 
scribed in Sec„ and treat the time dependent part V a (t) and 
15(f) as a perturbation. As a result, the effective self-energy 
E, which is due to the ac bias, would have two sources: the 
perturbation in leads Ejj, and the induced Coulomb interaction 
in molecular device U(/). Then, 



G\t,t') = G r ' ex (t 



J -o 



dh G r " eJC (f,fi)U(fi)G r (fi,0 



dh dt 2 G nex (t,h) 



G\t 2 ,t') 



where X5(t) is the response of the molecular device that is due 
to the Coulomb interaction when the time-dependent voltage 
is turned on. Here we have assumed an adiabatic response 
since most of time the variance of the applied electric field 
is much slower than the particles' intrinsic lifetime inside the 
scattering region. Then we have 15(f) = ±\50(-t) for down- 
ward case and upward case with U = H^f - Hj?. 

fdtidt 2 = lf dt\ f dt 2 + f dt\ f dt 2 \ 

%J \%J — oo %J — oo <J0 *J —oo i 



%(t,t') = V a (t,t')-T%\t-t') 
S£"(r-O = 8i(t)55 (t-O»««') 



Exact expression of kp(t, e) and F^f, e) 
Following the derivations in 

RefS we can get the exact 
expression for Ag(f , e) and Pp a (t, e) with the aid of the expres- 



sions £p — e + eVp and ep a — e + eVp - eV a : 

Ap (t, €) = G (e) + I — e 

x G r '°(E)[z(e (S )-Z(e)+P D G r ' v (e^] (29) 

F^(f,e) = J ^F(e)G«>\e)-Lf(E) + J g e « 

x \[z\e p ) - Z\e) + G a ' v (e fi )Pl] G aX) (E)Q D (E) 
+ [r( £ p a )G a - v (e fS )-Z*(e)G a -°(e)]i: a f(E)} (30) 

A%(t,e) = G^(*) + j g 

x G T ' v (E)[Z(e)-Z(e p ) + P u G r -\e)\ (31) 

F£(f,e) = J gz*(6^)G^)Sf(E) + Jg e -^ 

x {[z*(e) - T(e p ) + G a -°(e)Pl] G a ' v (E)Qu(E) 
+ e ieV °> [z*(e)G fl -°(e) - Z* (^ a )G a ' v (e p )\ 

(32) 



where 



P D = Z(ep)\] + Y J [2{ep)-Z(e ps )\[^\ep 6 )-i:f(E)\ 

s 

P v = -Z(e)U + ^[Z(e)-Z(Q)][Ef(e)-Ef(£-y a )] 

Qd(E) = J g [l-e^'- E »]Z(e')i:f(e') 

Qu(E) = J g [l Z(<)I£V) (33) 



with 



+M-1 



Z(e) = [*(£-£- *0 + )] 



(34) 



In the absence of the ac bias, the quantity A a is the Fourier 
transform of the retarded Green's function while the quantity 
Fp a is related to the Fourier transform of the advanced Green's 
function. They are all expressed in terms of the unperturbed 
Green's functions G r ^ a '°' v and self energy E°^ v which have 
been widely studied in molecular device using the NEGF- 
DFT formalism. G'^ a,0 ^ v and self energy E°^ y can be ex- 
pressed as 

G r/a,0/V {€) = [ £l _ H 0/V _ ^MO/^]" 1 

^•V) = [T° a -6Sy g L(e)[T ac -6Sy 
J7f(e) = [T° ca - <?S° ca ] gl a (e a ) [T° ac - e"S° ac ] 



where y 



r, a, <, e° 



e - eV a . Obviously, E„ (e) 



(e - eVa). In the wideband limit, Eq. (l29H321 > will reduce 
to the formula first derived by Jauho et al J25ll With A and 
F obtained we can, in principle, solve the AC current biased 
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by downwards or upwards pulse exactly. In practice, how- 
ever, its computational cost is expensive for a realistic molec- 
ular device. For example, to calculate J°'"(t), we have to do 
triple integrals over energy and repeat this procedure to collect 
data for all time sequence. In the numerical calculation espe- 
cially in ab-initio modeling, it is practically very difficult if 
not impossible to calculate the transient current for the com- 
plex structure in molecular devices. So approximation must 
be made so that Eq. (1291132b can be simplified. 



Approximate scheme of Ap(t, e) and Fp a (t, e) 

The approximate solution of Ap{t, e) and Fp a {t, e) in Eq. 
l32l have to satisfy the following requirements. First, it has to 
greatly reduce the calculational cost. Second, it has to keep es- 
sential physics of transient dynamics. Third, it must have the 
correct initial current at t — and approach the correct asymp- 
totic limit at t — > oo. The first goal is realized by eliminating 
double energy integral using a reasonable ansatz, with which 
the dynamical properties of molecular device is maintained. 

To find such an ansatz, we first assume that H a,0 {E) 
changes smoothly and slightly with E and is analytic in 
the upper half plane, so that the typical integral like 
j dedE zj^z^Qf ; S fl,0 (£') is roughly zero due to the different 
phase in e 1(e ~ £) '. Then the last term of F U ^ D and the second 
term of Q U ^ D disappear. Considering the following identity, 



duce the effective Green's function 



dE 

2^ -i{E -e + i0+) 
Or 1 rO 



r 1 f n CdE e' hT 

f-\f 

J-co * JO- 



dre i "'L a a (T) = 'L a a (e)-\ a a 



and defining A) = 2£(E) - A°, the first term of F u/D 

and Qu/d in Eqs.d33l> can be simplified, Fu/d now becomes 

F° * G a - (6)Sf( ej A) + j §-e-^< 

x [z^)-Z*(e) + G a ' v (ep)¥ f D ]G a - (E)i:f(E,A) (35) 

x [Z*(e) - T(e p ) + G a >°{e)Pl] G a ' v (E)2%°(E - eV a , A) 

(36) 



We note that, in the wide-band limit, Eq.( l35l36b is exact. With 
our approximation we have eliminated one of the energy in- 
tegrals in J'"", and A and F now have similar structures since 
F ~ A'V. 

With the approximation defined in Eq. (l35l36t . the current 
can be written in a compact form (see section C) if we intro- 



G r/a '°{E, e) 



G r,a ' V (E, e) 



-l 



ES 



(37) 



(38) 



In general we have to consider the overlap matrix S. How- 
ever, we should keep in mind that in the deriving of the 
time dependent current, we have to orthogonalize the basis 
set, which would lead to S = I. Here G r ' a (E, e) can be 
regarded as the Green's functions at energy E and constant 
parameter e for open system with the effective Hamiltonian 
Wj" = H c + X£(e). For a given H e// , Eqs. (l37l38b are equiva- 
lent to 

{ES - U r eff )G r - I (39) 

On the other hand, Green's function can be expanded in terms 
of the eigenfunctions of the corresponding Hamiltonian, 



G' = V *P"C„. 



(40) 



where H e// T" = E„(e)T„. Substituting Eq.(l40l) into Eq.(|39l, 
and using the general orthogonality relation 0" f S*P m = 
C m 6 nm [see Appendix J and the eigenvalue equation H e //T" = 
E n {ey¥", we have 

Z«F"0"t 
; (41) 
[E - E„(e)]0".tST" 

Obviously, this Green's function can be calculated by finding 
the residues Res,, = V^ 'ST" at various poles E = 

E n {e). 

Then, we replace Z(e)G r/a (E) in Eqs. d29|31l35l36l by 
Z(e)G r/a (E, e). Although G r,a {E, e) is different from initial 

Green's function G r,a {E) = [e - H c - ?. r l a {E)\\ this sub- 
stitution is reasonable since the major contribution of the in- 
tegration in Eqs.( |29ir32b comes from the pole e in Z(e) (see 
Eq.j34t). Similarly, considering the major contribution of 
the pole of Z(e), we replace Z(e)E fl -°(£) in Eqs. J29l3 1135136b 
by Z(e)E fl,0 (e). Since £(e) in G r (E, e) is independent of en- 
ergy E, we can perform contour integration over energy E in 
Eqs.d29b and OTI ) by closing a contour on lower half plane and 
perform the integration over energy E in Eqs.(|3"0Tl and ( f32b by 
closing a contour on upper half plane. Thus, energy integra- 
tion over E can be analytically performed. It should be noted 
that the self energy YI^ a is not independent of energy in con- 
trast to the wide-band limit, this energy dependence is on e 
but not on E. In this way, we can reduce the computational 
cost and keep the essential physics of the dynamics as we will 
show later. 



Approximate expression of A^f, e) and F fia (t, e) 

Now, considering the initial current and the asymptotic long 
time limit, we can write the approximate expression of A j g(r, e) 
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and Fftait, e) from Eqs. d29l31l35l36t : 

A D/U (t el - A D/U +A D/U 



(42) 



Fg,(f, e) = Aj^Ig ^, A) + A^'sf (e, A) (43) 
F£(f, e) - A^Ef (e, A) + A^X£° A) (44) 



with 



A° 2 = G r %) - J g e*^ [z(e)G- (£, e)] (46) 
A^ = J g e ^[Z(6)G^,£)(l + S^ (£))}47) 
A£ = G r -\ep) - J g [Z(^)G^(E, *)] (48) 

where 

a" = u + X[sfM-sfM] 

<5 
(5 

B D = -U + ^[sf(6)-Ef(e-^)] 

<5 

= -U + ^[Lf(e)-Ef(e)] (49) 

(5 

This is the second level of approximation. As we will see later 
that it is better than the first level approximation described be- 
low. Now we can make further approximation (the first level). 
To do this, we note that the Green's function G r can be ob- 
tained using the Dyson equation, 



where G r '"" is the Green's function of system denoted by 
H' of = H" + H', G r,ex is the unperturbed Green's function 
corresponding to H" that can be exactly solved, S is the ef- 
fective self energy describing H'. If we set BP and H'°' as 
zero biased open system and V a biased open system respec- 
tively, we have 

G r,M = G r,V {£) = G r,0 (e) + Qr,0 {e) -D - Q r,V ^ ^ 

Similarly, if we treat H ex and H'°' as V a biased open system 
and zero biased open system, respectively, we obtain another 
Dyson equation 



G r,tot = G r,u (e) = G r,v (£) + G'^e^C'V) (52) 

Similar to the derivation of the second level of approximation, 
we can also replace G"(e) by &"(£, e) in Eq. ( 151 1521 which 
leads to 

G nV (E, e) * G rX \E, e) [i + E D G r - v (e)] 

G r '°(E, e) * G'~- V (E, e) [i + S u G r -°(e)] (53) 



Then, Eqs.(l45ll and d47T > can be further approximated as 

A£ = J ^e^'[z( £ )G><\E,e)] (55) 

This is the first level of approximation. It is easy to confirm 
that when the self-energy is energy independent these two ap- 
proximations lead to exactly the same expression of transient 
current in the wide-band limit. In the next section we will nu- 
merically compare these two approximations with the exact 
solution. 

initial and asymptotic currents 

We now show that the currents calculated from 
Eqs. (123l28l42ll48l l and from Eg s . d2 3 128 1421144146 148 154|55l > 

satisfy the correct current limit at initial t — and asymptotic 
limit t — > c» times. Note that the initial current and asymp- 
totic currents can be calculated from a standard DC transport 
nonequilibrium Green's function analysis. It is expected that 
the asymptotic current for the downward pulse J%(t — » oo) 
and initial current for the upward pulse J^(t = 0) are zero 
since there is no bias in the system. Now we discuss the 
limiting cases for two versions of approximations developed 
in section IIIC. 

When t — 0, e'^~ E>t — 1, we can perform integration over 
energy E in Eqs. d45ll4"8T l by closing a contour at upper half 
plane, where only a single residual exists at an energy pole of 
Z. At t = 0, G r/a ' 0/v (E, e) = G r/a ' 0/v (e), therefore Eqs. d43R7l 
and Eqs.( l54l53T l are equivalent. Now we focus on the current 
obtained from Eqs. d23l28l42ll44l46l48l54l55b . After integrat- 
ing over e, the two terms in Eqs. fl46|48l l cancels to each other, 



(50) then from Eq. 



, A®' u (t = 0) becomes 



A£(f = 0) = G r ' v {ep,ep) = G r ' v {ep 
A%(t = 0) = G r ' (e,e) = G r -°(e) 



(56) 
(57) 



For Fp a , we can perform integration over energy E by closing 
a contour at lower half plane. Similarly, there also exists only 
a single residual on energy pole E z of Z* in the lower half 
plane, and 



F£(/ = 0) = G a - V {e p , epfE^iepa, A) = G a - V {e p )?% v {ep, A) 

(58) 

F*(t = 0) = G a - V {e, A) = G fl -°(e)£f (e, A) (59) 

Substituting Eq. d56ll59l ) into Eq. ( 1231281 ), and considering 



?.f(e) = (e,) 



:<,o/v 



(e) = G^V) 



ia,0/V 



E;- (£) = /(e)[Ef(6)-^°(e)] 

E|' v (e) = Ke - eV p ) [x a /(e) - U p v (e)\ (60) 
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where /(e) is Fermi distribution function, we have initial cur- 
rent at t — 

J° = 2eRe j ^ G r - V (e)5£ v (e) + G < - V (6)£j v (e)(61) 

J v a = 2 e Re J g G r ' (e)E<- (e) + G < -°(e)^°(£) (62) 

Eqs.doTTl and (l62l are the same as the formal DC current ex- 
pression in the case of nonzero bias and zero bias, respec- 
tively. jU(t = 0) in Eg. (1621 is exactly zero since the Fermi 
distribution in and G* are equal for a - L and a — R. 

When t — > oo, by virtue of the Riemann-Lebesgue 
lemma, [34] the Fourier integral over e vanishes, i.e., 
J ^ e _K *G r S7... equal to zero at t — > oo since there always 
exist poles in lower half plane. With this in mind, we have 



A£(r^oo,e) = G r '°(e) 



(63) 
(64) 

*8 V - ' — vj vcg; (65) 

F£(/ -» oo,e) = G°' y (^)Ef {epa, A) = G fl - V (6^ V (6 5 ,A) 

(66) 



F£(f-»°o,e) = G fl ' u (e)I£V,A) 
A^-»oo,e) = G r ' v (^) 



From Eq. d63ll66T ) and Eq. d23l28b . we have the asymptotic cur- 
rent 

] D a = 2 e Re J g G r ' (e)E<- (e) + G < -°(e)^°(e) (67) 
J v a = 2 e Re J g G r ' v (e)E<- v (e) + G<- V (6)S£ V (6) 

(68) 

It is easy to see, Eqs.doTb and d68l l are the formal DC current 
expression in the case of zero bias and nonzero bias, respec- 
tively, and J%(t — » oo) in Eq.d67|i is exactly zero. 



COMPARISON WITH THE EXACT RESULT IN QUANTUM 
DOT SYSTEM 



Now we consider a system composed of a single-level 
quantum dot connected to external leads with a Lorentzian 
linewidth. This system can be solved exactly to give a tran- 
sient current for pulse-like biasjfj. We can obtain transient 
current using three methods: (i) the exact current expressed 
by Eqs. d23|2"8l [29ll32l , (ii) the first level of approximation 
from Eqs. d23l28l42H44l46l48l54"l55| > and (iii) the second level 
of approximation from Eqs. d23l28l42ll48l >. We will compare 
the current obtained from these three methods. The system is 
described by the following simple Hamiltonian 

H = Y l ek a (t)4 a c ka + e d {t)d'd + Yj^A„ d + h - c -) < 69 ) 

k„ k„ 

where e d {f) = e" + U(t) and e ka (t) = e° + V a (t). Because 
the scattering region has only one state with energy level e", 



the Green's functions G(e) and self energy Z(e) thus become 
scalars instead of matrices. If we choose linewidth function 
r a (u>) = 2np a (u))\ti; a | 2 to be Lorentzian with the linewidth am- 
plitude r°, 



w 2 



then G 7 (e) and 27(e) can be expressed as 



G r/a '°(e) : 
G r/a ' v (e) 



e-e»-^Z'/-°(e) 

a 

e-eO-t^-^I'^e) 



;<,0/V 



;<,0/V 



U) 



. a 

sf- (6) = J ^r o (c)/( £ -u±i0 + ) 

E<-°(e) = /(e) [Ef (6) - Sf (e)] 
2£ v (e) = /(e - eV a ) [l£ y (e) - (6)] 

Using the theorem of residual, we can analytically perform 
integral in Ap and Fp a for either exact formula or two approx- 
imate formulas. In the calculation, we set F = F" + T° R as the 
energy unit, and set T° L = T° = 0.5. 

We first consider the transient current induced by opposite 
voltage Vi(t) = -Vfi(f)- In this case, the equilibrium coulomb 
potential in quantum dot U 0,v = 0, and the time dependent 
perturbation coming from coulomb response U (f) is assumed 
to be zero. It is a reasonable assumption since the coulomb po- 
tential in scattering region is canceled by the opposite voltage 
in left and right lead. In Fig. 1, we plot two approximated tran- 
sient currents and exact transient current in downward [panel 
(a), (b), (c)] and upward [panel (d), (e), (f)] case vs time for 
different bandwidth W. We find that for all bandwidth W, the 
approximated current and exact current have the same dynam- 
ical behaviors. Fig. 2 gives direct comparison where we merge 
panels (a), (b) and (c) in Fig.l as panel (a) in Fig. 2, and merge 
panels (d), (e) and (f) in Fig. 1 as panel (b) in Fig. 2. We can see 
that for the downward pulse [panel (a)], transient current us- 
ing three formulas are almost indistinguishable. This means 
that in the opposite voltage, our approximation, the first ap- 
proximation [Eqs. (l46l48l54l55H and the second approxima- 
tion rEqs. (!45H48l >l are all very good for studying transient dy- 
namics. For the upward pulse, although the approximations 
are not as good as in downward case, the currents calculated 
from approximate scheme are still in good agreement with the 
exact solution especially for the second approximation. Hence 
we may conclude that the two approximations are all reason- 
able in the opposite voltage VzXO = -V«(0- They can be used 
to study transient dynamics in the real molecular device to 
speed up the calculation. 
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FIG. 1: (Color online) Time dependent current J{t) corresponding 
to an opposite downward pulse or upward pulse in three versions: 
the exact solution and two approximations. The different black lines 
are for different bandwidth W. The red line is for W = °°, i.e., the 
wide-band limit. The current is in the unit of eT, and the time is in 
the unit of 2n/T, eV L = -eV R = 5. 



Next, we focus on the asymmetric voltage, i.e., Vl(0 + 
VR(t). In this case, the equilibrium coulomb potential in quan- 
tum dot U°t y , and the time dependent perturbation coming 
from coulomb response U(t) can't be canceled by the voltage 
in left and right lead. In principle, perturbation U (t) should be 
calculated by solving time dependent Schrodinger equation, it 
will be very difficult and computational demanding therefore 
can't be implemented in real molecular device. As an alter- 
native scheme, we have set U(t) - [eV L (t)Y° L + eV L (t)r° L ]/T. 
For the single level quantum dot system, this is exact because 
the central scattering region now is expressed in a scalar in- 
stead of matrices, which leads to the same transient current 
for the opposite voltage V/,(/) = -Vr(£) and asymmetric volt- 
age V L (t) = V(t), V R (t) = or V L {t) = 0, V R (t) = -V(t) in the 
exact solution. 

For the first approximation the poles in time dependent term 
e At-E)t are different from that in the second level approxima- 
tion, i.e., the poles of G''° in Eq|45]and G r ' v in EqgT] are re- 
placed by the poles of G r,v in Eq|54] and G' ° in Eq|55] re- 
spectively. Because of this, the time evolution process are not 
as accurate in the first approximation, especially for the large 
V a . So, for the asymmetric voltage, the second approxima- 
tion is better. In Fig. 3 and Fig. 4, we compare the transient 
current obtained from the second approximation [panel (b-d)] 
for opposite or asymmetric voltage with the exact transient 
current [panel (a)] in response to the downward pulse and up- 
ward pulse, respectively. We find that all transient currents 
from the second approximation in Fig. 3 and Fig.4 [panel (b)] 
are very close to the exact result [panel (a)]. Moreover, in 
Fig. 3 and Fig.4, the approximate transient current in panel (b), 
(c), (d) have almost the same behavior. It is safe to say that 
our approximations have kept essential physics of dynamical 
transport properties. 



FIG. 2: (Color online) Merged version of Fig.l for W = 1, 2, 5 and 
20. Panel (a) corresponding to the downward pulse current comes 
from panel (a), (b) and (c) in Fig.l, panel (b) corresponding to up- 
ward pulse current comes from panel (d), (e) and (f). Along the black 
arrow, the bandwidth are W = 1, 2, 5 and 20, respectively. 



SEVERAL EXAMPLES FOR REAL MOLECULAR 
DEVICES 

In this section, we implement our approximate formula in 
two representative molecular devices including a short carbon 
chain coupled to aluminum leads and a Ceo molecule coupled 
to aluminum leads. These systems were chosen because they 
are typical in first-principles calculation and their practical 
importance to nano-electronics. In Fig. 5(a) and Fig. 5(b), we 
show the structure of AI-C5-AI and AI-C60-AI, respectively, 
where Al leads are along (100) direction, one unit cell of Al 
lead consists of 9 Al atoms and total 40 atoms were included 
in the simulation box. For the A1-C5-A1 device, the nearest 
distance between Al leads and the carbon chain is 3.781 a.u. 
and the distance of C-C bond is 2.5 a.u.(l a.u. =0.529 A). In 
A1-C60-A1 device, the distance between the Al atom and the 
nearest C atom equal to 3.625 a.u.. 

To calculate the dynamic response of molecular devices, 
we have used the first-principles quantum transport package 
MATDC AL. J35I Considering the complicated coulomb re- 
sponse in scattering region, we set V[£t) = -Vr(0- in this 
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FIG. 5: (Color online) Panel (a): Structure of A1-C5-A1. Panel(b): 
structure of A1-C 60 -A1. 



FIG. 3: (Color online) Panel (a): exact time dependent current J(t) 
corresponding to downward pulse for dV = V L — V R = 5. Panel 
(b-d) are corresponding to the second approximate transient current 
corresponding to downward pulse for opposite voltage V L = - V R = 
2.5, asymmetric voltage V L = 5, V R = and V L = 0, Vr = -5, 
respectively. The different black lines are for different bandwidths 
W. The red line is wide-band limit for W = oo. 



FIG. 4: (Color online) Same to Fig(3] transient current corresponding 
to upward pulse vs time are plotted. 



case, the first approximation is simple but as good as the 
second one. So, in the following, the first approximate for- 
mula [Eqs.( l23l28l42H44l46l48l54T55l ll is used. In principle, 
the calculation involves the following steps: (1) calculate the 
device Hamiltonian including central scattering Hamiltonian 
and lead Hamiltonian using NEGF-DFT package to get two 
potential landscapes U at zero bias and U v at V a bias, re- 
spectively. They are originally expressed in a nonorthogo- 
nal fireball basis. (2) orthogonalize the nonorthogonal de- 
vice Hamiltonian using the approachj3^| introduced in Ap- 
pendix„so that they are finally expressed in an orthogonal ba- 
sis. (3) with the orthogonal lead Hamiltonian H a , one calcu- 
lates zero biased self energy "L'J"' and V a biased self energy 
~L r J, a,v from Eqs. (ll3ll4l i using the transfer matrix method.! 31 ] 



(4) with orthogonalized central scattering Hamiltonian H[! and 
and self energy and H'J a ' v obtained from two poten- 
tial landscapes U° and U v , one solves the effective Green's 
function Q r l a '°l v using Eqs. (|37|38t by calculating its poles 
and residuals from Eq.flTb. Step (l)-(4) are time independent 
processes and easy to perform. (5) calculate time dependent 



Then hp and Fp a can be calculated from Eqs.( 142H4"4l ). (6) 



quantities Ag[ u and A^' 2 U from Eqs. (l54l55l and Eqs. (l46l48l >. 



integrate over e and obtain the final AC current J D ^ u {t) = 
[J° /U (t) - J% /U (t)]/2 from Eqs. (l23l28l . 

First we study the AI-C5-AI structure. In Fig. 6, we plot 
the transient current J(t) corresponding to the upward pulse 
[panel (a) and (b)] and the downward pulse [panel (c) and (d)] 
for different external voltages Vr = -Vl = 0.001a. u. [panel 
(a) and (c)] and V R = -V L = 0.01a.m. [panel (b) and (d)] in Al- 
C5-AI structure. Following observations are in order: (1) as 
we have discussed in Sec„ for all bias voltages V a the transient 
currents indeed reach the correct limits at t — and t — > 00. 
For the upward pulse, J(t — 0) = and J(t — > 00) = J dc 
while for the downward pulse we have J(t = 0) = Jdc and 
J(t — > 00) = 0. (2) for both upward pulse (turn-on voltage) 
and downward pulse (turn-off voltage), once the bias volt- 
age is switched, currents oscillate rapidly in the first a few or 
tens fs and then gradually approach to the steady-state values, 
i.e., Jdc for turn-on voltage and zero for turn-off voltage. The 
larger the voltage V a , the more rapid the current oscillates. (3) 
concerning the long time behavior, the time dependent current 
oscillates with a frequency proportional to | V a \ . 112.211 This is be- 
cause the time dependent term e' (e ~ £) ' in Eqs. d46|48l54|55l ) are 
V„ dependent. For the upward pulse, oc e' v °', which 

directly leads to the oscillating frequency proportional to \ V a \. 
For the downward pulse, although g'( £ ~ £)/ is V a independent, 
in the energy integral on E, the pole E„ of G' (E, e) are de- 



termined by the self energy . Since depends on V a , 
this leads to V a dependent oscillating frequency. In addition, 
we notice that although the properties of dc conductance of 
short carbon chains are different for the chains with odd and 



even number atoms 13711 due to the completely different elec- 
tronic structure near Fermi level, the ac signals are similar 
(see Ref I22I where AI-C4-AI structure was analyzed). This in- 
dicates that in AC transport, all states with energy from -00 
to the Fermi energy are contributing, which is very different 
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FIG. 6: (Color online) Time dependent current J(t) corresponding to 
the upward pulse [panel (a) and (b)] and the downward pulse [panel 
(c) and (d)] for different external voltages V a in A1-C 5 -A1 device. The 
inset of panel (a) shows the long time behavior of the time-dependent 
current. The red (gray in print) dashed lines in panels indicate asymp- 
totic current J(t — » 00) which the DC current biased by Vl/r labeled 
in corresponding panels for the upward pulse, and arrive at zero for 
the downward case. 

from dc case where only the states near Fermi level contribute 
to transport processes. 

Next, we study the second sample: the AI-C60-AI structure. 
In Fig. 7, the transient current J{f) of the structure correspond- 
ing to an upward pulse [panel (a) and (b)] and a downward 
pulse [panel (c) and (d)] for different external voltages Vr = 
-V L = 0.001a.m. [panel (a) and (c)] and V R = -V L = 0.01a.m. 
[panel (b) and (d)] are plotted. Similar to the AI-C5-AI struc- 
ture, correct initial current J(t = 0) and asymptotic current 
J(t — > 00) are also obtained in AI-C60-AI structure. In addi- 
tion, there are also rapidly oscillations at short times after the 
switch although the oscillation is not as rapid as that in the Al- 
C5-AI structure. Furthermore, similar to AI-C5-AI structure, 
in gradually reaching the steady-state values, the current os- 
cillates with a frequency proportional to \ V a \ but its decay rate 
is much slower than that in AI-C5-AI structure. It indicates 
that there are much more quasi-resonant state that contribute 
to the transient current in Al-Ceo-Al structure which is reason- 
able considering the complex electronic structure of isolated 
Ceo- In the following, we will analyze in detail how the cur- 
rent decays for the AI-C60-AI structure. 

Physically, decay time of current corresponds to the width 
of the quasi-bound state. In molecular devices, because the 
linewidth function T(e) are complex and energy dependent 
matrix, we can't extract characteristic time scale directly from 
1 /r. As such, the transmission coefficient T(e) is needed to 
understand the resonant state and corresponding characteris- 
tic time scale. In Fig. 8(a), we plot transmission coefficient 
T(e) in the energy range from the energy band bottom to the 
Fermi energy for A1-C60-A1 structure at zero bias. Here, the 
sharp peaks [some of them, see red crossed signed peaks in 
Fig. 8(a)] correspond to resonant states with large lifetimes. 



FIG. 7: (Color online) Time dependent current J(t) corresponding to 
the upward pulse [panel (a) and (b)] and the downward pulse [panel 
(c) and (d)] in AI-C60-AI device for different V a . In panel (a) and (c), 
V R = -V L = O.OOla.H.. In panel (b) and (d), V R = -V L = O.Ola.u.. 
Same to Fig. 6, the red (gray in print) dashed lines in panels indicate 
asymptotic current J{t — > 00). The long time AC current or detailed 
short time AC current are shown in inset of panels. 



At a particular resonant state, the incoming electron can dwell 
for a long time, which contributes to a much more slowly de- 
caying current than other non-resonant states. In Fig. 8(b), (c) 
and (d), we amplify the first, second and forth labeled quasi- 
resonant transmission, respectively, where the peaks' width 
^peak ~ 10 -5 a.w. are indicated, corresponding to a decay time 
t ~ 2400/ s from the expression T pea kt =1. In Fig. 8(e)- 
(g), corresponding to different e where the resonant peaks in 
Fig.8(b)-(d) are located, we plot long time behavior of current 
element 7z.(e). Here J^e) is the time dependent current for 
each energy e, the integration over which gives the final cur- 
rent J a (t). We can see that for each resonant state the current 
Jl(£) keeps oscillating in a long time comparable to the decay 
time t ~ 2400/s. Furthermore the intensity of the oscilla- 
tion AJ ~ 0.2/jA is not very small comparing to the DC signal 
Jdc = 5.lfiA. 

After integration over energy, these slowly decaying cur- 
rents Jl(£) due to the resonant states may cancel to each 
other partially due to the difference in their phases. However, 
we should keep in mind that it is these resonant peaks that 
may give rise to convergence problem. Hence in the calcula- 
tion, we should first scan the equilibrium and non-equilibrium 
transmission coefficient (100,000 energy points for example) 
to resolve sharp resonant peaks in the whole energy range 
from minimum energy to Fermi energy. Then, for each sharp 
resonant peak, enough (100 for example) energy points should 
be chosen to converge the integration of the current Jl(s) over 
e, i.e., J deJ(e). For the non-resonant state, i.e., the smoothly 
changed region in T(e), the current J(e) are integrated using 
less energy points. 

As we have discussed that the resonant states are important 
for the transient current and they must be carefully treated 
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FIG. 8: (Color online) Panel (a): transmission coefficient T(e) in the 
energy range from the energy band bottom to the Fermi energy. In 
the whole energy range, there are some resonant states corresponding 
to the very sharp transmission coefficient T(e), as we have indicated 
(see red cross) and labeled (by 1, 2, 3 and 4) in panel (a), some of 
them contribute to the current at long time. We amplify the first, 
second and forth labeled resonant transmission in panel (b), (c) and 
(d), respectively. In panel (e)-(g), we plot the long time behavior of 
current Jl(c) at a fixed e for the first, second and forth resonant states. 
The external voltage V a = 0.001a. z<.. 



in calculation. However, in the calculation of the effective 
Green's function G r ^ a '°^ v , a small imaginary part that is usu- 
ally added to the real energy e — > e + it] to help resolving 
the retarded or advanced self-energies. This in turn introduces 
pseudo resonant states. In order to eliminate the pseudo reso- 
nant state in effective Green's function G'~ /a ' 0,v [Eqs. d37|38"l )|, 
one has to calculate the self-energy by setting 77 = and re- 
solve the retarded or advanced self-energies with the aid of the 
group velocity v* — (dE(k)/dk).]2&\ 



CONCLUSION 

By orthogonalizing the Hamiltonian expressed in the 
nonorthogonal basis and considering the singularity of self- 
energy Y, r l a (t, t') at t — f, we have generalized the solution [ 
developed in RefH of the transient current driven either by 
a downward step voltage pulse or by a upward step pulse. 
This generalized result can be applied to both the quantum 
dot model and real molecular device. Based on the exact solu- 
tion given in Ref.6, we derived two approximate formulas that 
are suitable for numerical calculation of the transient current 
for molecular devices. We have tested our approximate for- 
mula in a quantum dot system where exact numerical solution 
exists. For the quantum dot system, we chose a Lorentzian 
linewidth (beyond wideband limit) and compared the time- 
dependent current calculated using both exact formula and our 
approximate formula. We found that for the opposite voltage 
Vl(/) = — Vr(0, the results obtained from the exact formalism 
and two approximate scheme agree very well with each other 



especially in the downward pulse case. For the nonsymmetric 
voltage V L (t) = V(t), V R (t) = or V L (t) = 0, V R (t) = -V(f), 
the second approximation is better. This shows that our ap- 
proximate formulas captured the essential physics of the tran- 
sient current. In addition, it gives the correct initial current 
at t = and correct asymptotic current at t — > 00. Since 
we have reduced the calculation from triple integral to single 
integral over the energy, the approximated approach reduces 
the computational cost drastically and it can be easily imple- 
mented in first principles calculation for molecular devices. 
To demonstrate this, we calculated the transient current us- 
ing the first approximated scheme with an opposite voltage 
V/XO = -Vu(t) for two molecular structures: AI-C5-AI and 
Al-Ceo-Al. Different from the quantum dot system, because 
of the complex electronic structure in molecular devices, tran- 
sient currents oscillate rapidly in the first a few or tens fs as 
the bias voltage is switched, then gradually approach to the 
steady-state values. Furthermore, due to the resonant state in 
molecular devices, transient currents have a very long decay 
time t. 



orthogonality relation for the nonorthogonal basis 

For a system described by H, the time independent eigen- 
value equation is written as: 



H\n) = E n \n) 



(70) 



the eigenvectors \n) form an orthogonal complete basis set. 
However, in many systems such as a molecular device con- 
nected to external leads, the basis set constructed by eigen 
vectors is not convenient. We usually expand the eigen vector 
\n) in other basis which is non-orthogonal complete set (or 
nearly complete). 

|n>^|/i><//|n> (71) 
the eigenvalue equation now becomes 

Y^vm^W) = £„2<v|iU><iu|n) 

n ft 

Y^^l = E„]TS VII V; (72) 

H ft 

where S v/J = (v|/i). In the matrix form, we have HT" = 
EnSV. If we use the self-energy to replace the effect of 
leads the effective Hamiltonian for the open system becomes 
H = Ho+57. Since the effective Hamiltonian is not Hermitian, 
we can define the adjoint operator H 1 = H = Ho + 2 a and cor- 
responding eigen-equation becomes H + |0„) = E*S\(f>„). Then 



0" ,+ H¥" = £ n O" ST", 
«P". t H + 0"' = ElW nf S f <S>" 



(73) 
(74) 
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Taking hermitian conjugate of Eq.d74l. 

cl) m HT" = £ m O" vl "ST" (75) 

From (T73T > and (fTBT l. we have 

0-.t SV « = Cm4m (76) 

For the normalized wave function \\j/ n ) and \<p„), 

OSI* = I (77) 

It is the usual orthogonality relation for eigenvectors ex- 
pressed in a nonorthogonal basis set. For an hermitian Hamil- 
tonian H = H 1 , \ifr n ) = \<p n ), we have 

*FS*F = I. 



Orthogonalize Hamiltonian expressed in nonorthogonal basis 

In this appendix, we will show how to construct a new or- 
thogonal basis from the atomic real-space nonorthogonal ba- 
sis. We will transform the original Hamiltonian H which is 
expressed in the nonorthogonal basis into Hamiltonian H ex- 
pressed in the new orthogonal basis. Of course, instead of S, 
the overlap matrix in the new basis will be I. 

Denoting nonorthogonal basis \fi) and orthogonal basis \ j), 
they are related by unitary transform operator H 



Up = <M 



(78) 



where we have used the completeness of orthogonal basis \ 
Using the orthogonality = 5,j 

YjAn)<4Av){v\j) = Yj<U ifl S tlv Ul j = 5 ij 

[IV fiV 

where we have used the completeness of nonorthogonal basis. 
In the matrix form, KSK^ = I. We can formally define 

Then new Hamiltonian H expressed in basis \i) can be ex- 
pressed as: 

H,7 = (W\j) 



(79) 



In the matrix form, H = S 5 H |S 2 J ' . 

We now discuss how to find the matrix S _I . Without loss 
generality, we assume the real overlap matrix S satisfies eigen 



function SV = Vdiag(/li, A„) with the eigenvalues A 1, A„ 
and eigenvectors V = [vi, v„]. Since S is real and symmet- 
ric, the eigenvectors are real and orthogonal, and it thus holds 
that V' V = VV 1 = I. Then 

S = Vdiag(^,...,/l„)V t 

= Vdiag( y[A x , V^)V" , "Vdiag( y/M, V^)V" f 



It follows that 



= Vdiag( V^)V"'" 



(80) 



From S * = I and Eq.(f80b. we have 

S^Vdiag( V^i, V^)V f = I 

S-Wdiag(V*T,..., V^^VdiagC-ip. 

V/ii 



)V f 



= S" 



: Vdiag( 



1 



1 



)V" f 



(81) 



In general, the dimension of matrix S is infinity, we can't 
calculate its eigenvalue A/ and eigenvector v, by diagonaliz- 
ing S. However, in the tight-binding representation, the state 
fi and v hardly overlap when their separation is large enough 
in real space, i.e., S^v « for most of off-diagonal elements. 
Considering the periodic properties in semi-infinite leads, we 
can select a block matrix which is large enough to include all 
the overlap between leads and central molecular regions. For 
the non-orthogonal basis including several unit cell of atomic 
leads as a buffer layer into the central scattering region is 
enough to get a good screening for dc transport calculation. In 
transforming the Hamiltonian to the orthogonal basis needed 
for ac transport calculation, however, it turns out that we have 
to include at least 10 unit cells of atomic leads into the central 
scattering region. Partly because the overlap of orthogonal ba- 
sis has longer range than that of non-orthogonal basis. With 
this large simulation box (finite dimension), we can calculate 
the overlap matrix S2 therefore transform H into H. The ac- 
curacy of transformed Hamiltonian H should be examined by 
comparing dc conductance obtained from the original Hamil- 
tonian H and the transformed Hamiltonian H. 
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